function [rej_rate,p_adj] = rej_rate_CRS_mean(Sigma_hat,membership,alpha_sig,alt)

rng(7)
n = size(Sigma_hat,1);
clusters = unique(membership);
G = numel(clusters);
nSim = 10000; 

% only generates group mean
L = repmat(membership',G,1)==repmat(clusters,1,n);
N_g_vec = sum(L,2);
L_weighted = L./repmat(N_g_vec,1,n);
b_mat = alt+(randn(nSim,G)*chol(L_weighted*Sigma_hat*L_weighted'))';
b_fm_vec = mean(b_mat); % FM estimates for all replications

nSim_crs = 1000; % number of simulated transformations
transformation_mat = (((rand(nSim_crs,G)>.5)-.5)*2)/G;

b_distri = transformation_mat*b_mat;

rej_rate = mean(abs(b_fm_vec)>quantile(abs(b_distri),1-alpha_sig)+eps);

% adjusted p-value threshold such that rejection rate is alpha_sig
pval_distr = mean(repmat(abs(b_fm_vec),nSim_crs,1)<=abs(b_distri));
p_adj = quantile(pval_distr,alpha_sig);



